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Nanoelectromechanical systems are characterized by an intimate connection between electronic 
and mechanical degrees of freedom. Due to the nanoscopic scale, current flowing through the 
system noticeably impacts the vibrational dynamics of the device, complementing the effect of 
the vibrational modes on the electronic dynamics. We employ the scattering matrix approach to 
quantum transport to develop a unified theory of nanoelectromechanical systems out of equilibrium. 
For a slow mechanical mode, the current can be obtained from the Landauer-Biittiker formula in 
the strictly adiabatic limit. The leading correction to the adiabatic limit reduces to Brouwer's 
formula for the current of a quantum pump in the absence of the bias voltage. The principal result 
of the present paper are scattering matrix expressions for the current-induced forces acting on the 
mechanical degrees of freedom. These forces control the Langevin dynamics of the mechanical 
modes. Specifically, we derive expressions for the (typically nonconservative) mean force, for the 
(possibly negative) damping force, an effective "Lorentz" force which exists even for time reversal 
invariant systems, and the fiuctuating Langevin force originating from Nyquist and shot noise of the 
current flow. We apply our general formalism to several simple models which illustrate the peculiar 
nature of the current-induced forces. Specifically, we find that in out of equilibrium situations the 
current induced forces can destabilize the mechanical vibrations and cause limit-cycle dynamics. 

I. INTRODUCTION 

Scattering theory has proved a highly successful method for treating coherent transport in mesoscopic systems [T]. 
Part of its appeal is rooted in its conceptual simplicity: transport through a mesoscopic object can be described 
in terms of transmission and reflection of electronic waves which are scattered by a potential. This approach was 
introduced by Landauer |2l3j and generalized by Biittiker et al. 0] and leads to their well-known formula for the 
conductance of multi-terminal mesoscopic conductors. For time-dependent phenomena, scattering matrix expressions 
have been obtained for quantum pumping |5|6j . a process by which a direct current is generated through temporal 
variations of relevant parameters of the system, such as a gate voltage or a magnetic field. The case of pumping in 
an out-of-equilibrium, biased system has remained largely unexplored so far |7I8I . 

The purpose of the present paper is to further develop the scattering matrix approach into a simple, unifying 
formalism to treat nanoelectromechanical systems (NEMS). The coupling between mechanical and electronic degrees 
of freedom is the defining characteristic of NEMS [9,10^ , such as suspended quantum dots pdj , carbon nanotubes or 
graphene sheets |12ll3j . one-dimensional wires [T3], and molecular junctions |15I16| . For these systems, a transport 
current can excite mechanical modes, and vice versa, the mechanical motion affects the transport current. The reduced 
size and high sensitivity of the resulting devices make them attractive for applications such as sensors of mass or 
charge, nanoscale motors, or switches [T7]. On a more fundamental level, the capability of cooling the system via 
back-action allows one to study quantum phenomena at the mesoscopic level, eventually reaching the quantum limit 
of measurement |18ll9j . 

All of these applications require an understanding of the mechanical forces that act on the nanoelectromechanical 
system in the presence of a transport current. These are referred to as current-induced forces, and have been observed 
in seminal experiments |20l21j . Recently we have shown that it is possible to fully express the current-induced forces 
in terms of a scattering matrix formalism, for arbitrary (albeit adiabatic) out of equilibrium situations |22l. thus 
providing the tools for a systematic approach to study the interplay between electronic and mechanical degrees of 
freedom in NEMS. 

In the context of NEMS, two well defined limits can be identified for which electronic and mechanical time scales 
decouple, and which give rise to different experimental phenomena. On one side, when the electronic time scales are 
slow compared with the mechanical vibrations, drastic consequences can be observed for the electronic transport, such 
as side bands due to phonon assisted tunneling |23l24j or the Frank-Condon blockade effect, a phononic analog of 
the Coulomb blockade in quantum dots In the opposite regime, electrons tunnel through the nanostructure 

rapidly, observing a quasistatic configuration of the vibrational modes, but affecting their dynamics profoundly at the 
same time [TBHH]. It is on this regime that our present work focuses. We treat the vibrational degrees of freedom as 
classical entities embedded in an electronic environment: pictorially, many electrons pass through the nanostructure 
during one vibrational period, impinging randomly on the modes. In this limit, it is natural to assume that the 
dynamics of the vibrational modes, represented by collective coordinates X^, will be governed by a set of coupled 



2 



Langevin equations 



dlJ 

M,X, + — =F,-Y, l^-'X.' + ■ (1) 



Here we have grouped the purely elastic contribution on the left hand side (LHS) of Eq. ([T]), M^, being the effective 
mass of mode v and ?7(X) an elastic potential. On the right hand side (RHS) we collected the current-induced forces: 
the mean force F^^, a term proportional to the velocity of the modes —^jji^vv'Xui, and the Langevin fluctuating 
forces . The main result of our work are expressions for the current- induced forces in terms of the scattering matrix 



and its parametric derivatives. These are given by Eq. (|39[) for the mean force F^{X), Eq. (42) for the correlator 



Diyi,'(X) of the stochastic force and Eqs. (47) and (|500or the two kinds of forces (dissipative friction force and 
effective "Lorentz" force, as we discuss below) encoded by the matrix 7^y'(X). 

Theoretically, these forces have been studied previously within different formalisms. The case of one electronic level 
coupled to one vibrational mode has been studied with a Green's function approach in Refs. 28 29J, where the authors 
showed that the current-induced forces can lead to a bistable effective potential and consequently to switching. In 
Ref. |30j , the authors studied the case of multiple vibrational modes within a linear approximation, finding a Lorentz- 
like current- induced force arising from the electronic Berry phase [3T]. In simple situations, the current-induced 
forces have been also studied within a scattering matrix approach in the context of quantum measurement backaction 
[32] (see also [33j), momentum transfer statistics [34], and of magnetic systems to describe Gilbert damping j35| . 
Current induced forces have been shown to be of relevance near mechanical instabilities [55H55] and to drive NEMS 
into instabilities and strong non-linear behavior |39H41j. Our formalism allows us to retain the nonlinearities of the 
problem, which is essential for even a qualitative description of the dynamics, while turning the problem of calculating 
the current-induced forces into a scattering problem for which standard techniques can be applied. 

In what follows we develop these ideas in detail, giving a thorough derivation of the expressions in terms of the 
scattering matrix for the current-induced forces found in Ref. |22| . and include several applications to specific systems. 
Moreover, we extend the theoretical results of Ref. [52] in two ways. We treat a general coupling between the collective 
modes X^, and the electrons, generalizing the linear coupling expressions obtained previously. We also allow for an 
arbitrary energy dependence in the hybridization between the leads and the quantum dot, allowing more flexibility 
for modeling real systems. In Section [IT] we introduce the theoretical model, and derive the equations of motion of the 
mechanical degrees of freedom starting from a microscopic Hamiltonian. We show how the Langevin equation, Eq. 
([T]), emerges naturally from a microscopic model when employing the non-equilibrium Born Oppenheimer (NEBO) 
approximation, appropriate for the limit of slow vibrational dynamics, and derive the current induced forces in terms 
of the microscopic parameters. In Section |III| we show that the current induced forces can be written in terms of 
parametric derivatives of the scattering matrix (S-matrix) of the system, and state general properties that can be 
derived from S-matrix symmetry considerations. In Section [IV| we complete the discussion of nanoelectromechanical 
systems in terms of scattering matrices by providing a corresponding expression for the charge current. In Section [V] 
we apply our formalism to simple models of increasing complexity, namely a single resonant level, a two-level model, 
and a two- level/two- mode model. We conclude in Section |Vl] For better readability, we have relegated part of some 
lengthy calculations to the Supplementary Material, together with a list of useful relations that are used throughout 
the main text. 

II. MICROSCOPIC DERIVATION OF THE LANGEVIN EQUATION 

A. Model 

We model the system as a mesoscopic quantum dot connected to multiple leads and coupled to vibrational degrees 
of freedom. Throughout this work we consider non- interacting electrons and we set h = 1. The Hamiltonian for the 
full system reads 

H = Hd + Hx+Hl + Ht, (2) 

where the different terms are introduced in the following. 

We describe the quantum dot by AI electronic levels coupled to N slow collective degrees of freedom X = 
{Xi, . . . , Xj^). This is contained in the dot's Hamiltonian 



/io(X) dm> (3) 



mm' 
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which describes the electronic levels of the dot and their dependence on the collective modes' coordinates X^, [v = 
1, . . . , N) by the hermitian M x M matrix /io(X). The operator (d) creates (annihilates) an electron in the dot and 
the indices m, m' (= 1, . . . , M) label the electronic levels. Note that here we generalize our previous results obtained 
for a linear coupling in X |22j . and allow ho to be a general function of X. Our analysis is valid for any coupling 
strength. The free evolution of the 'mechanical' degrees of freedom of the dot is described by the Hamiltonian 



p2 

2M„ 



C/(X) 



The leads act as electronic reservoirs kept at fixed chemical potentials /Iq and are described by 



(4) 



(5) 



where we represent the electrons in the leads by the creation (annihilation) operators (c). The leads' electrons 



obey the Fermi-Dirac distribution fa (e) = [l 



,(£-M=)/fcTl 



The leads are labeled by a = 1, . . . , L, each containing 



channels n — 1, . . . , Na- We combine 77 — (a, n) into a general 'lead' index, rj — 1, . . . , Nq with A^o = J2a 
Finally, the Hamiltonian H-r represents the tunneling between the leads and the levels in the dot. 



Ht = ^(c^W^mdm + h-c.) 



(6) 



B. Non-equilibrium Born-Oppenheimer approximation 

We use as a starting point the Heisenberg equations of motion for the mechanical modes which can be cast as 

dU 



where we have introduced the X-dependent matrices 

dho 



Y,di\A.{±)] d„', (7) 



A 



The RHS of ([t]) contains the current-induced forces, expressed through the electronic operators d of the quantum 
dot. We now proceed to calculate these forces within a non-equilibrium Born-Oppenheimer (NEBO) approximation, 
in which the dynamics of the collective modes is assumed slow. In this limit, we can treat the mechanical degrees of 
freedom as classical, acting as a slow classical field on the fast electronic dynamics. 

The NEBO approximation consists of averaging the RHS of Eq. ([t]) over times long compared to the electronic 
time scale, but short in terms of the oscillator dynamics. In this approximation, the force operator is represented by 
its (average) expectation value ((i^A(i)x(t), evaluated for a given trajectory X(i) of the mechanical degrees of freedom, 
plus fluctuations containing both Johnson-Nyquist and shot noise. These fluctuations give rise to a Langevin force 
Hence Eq. Q becomes 

M,X, + ^ - tr[zA,5< (t, t)] + , (9) 
where the trace "tr" is taken over the dot levels, and we have introduced the lesser Green's function 

QL'M^MAt')dn{t))Mt)- (10) 

The variance of the stochastic force is governed by the symmetrized fluctuations of the operator S Kd. Given that 
the electronic fluctuations happen on short time scales, is locally correlated in time, 

{Ut)i.'(t'))^D,,,{iq5{t~t'). (11) 

(An alternative but equivalent derivation, is based on a saddle point approximation for the Keldysh action, see e.g. 
Ref. 142] ). Since we are dealing with non-interacting electrons, D(X.) can be expressed in terms of single particle 
Green's functions using Wick's theorem. This readily yields 

iUm.'it')) ^ tT{A,g>{t,t')A,,g<it' ,t)}, , (12) 
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where 

C^'(i,^') = -*(rf™WrfL(i')>x(t) (13) 

is the greater Green's function. These expressions for the current-induced forces show that we need to evaluate the 
electronic Green's function for a given classical trajectory X(i). In doing so, we can exploit that the mechanical 
degrees of freedom are assumed to be slow compared to the electrons. Thus, we can approximate the Green's function 
by its solution to first order in the velocities X(i). We now proceed with this derivation, starting with the Dyson 
equation for the retarded Green's function 

5Ln'M = -^Oit-t'){{d„,{t),dl,{t')})^^t). (14) 

Here {., .} indicates the anti-commutator. We note that since we consider non-interacting electrons, we can restore the 
lesser and greater Green's functions (or the advanced Green's function Q^) at the end of the calculation by standard 
manipulations. 

The hybridization with the leads is taken into account through the self-energy [43j 

j:%)^-zJ2r^ie), (15) 

a 

which is given in terms of the width functions 

r^{e)^nW^{e)Il„W{e). (16) 

Here we have defined Hq, as a projection operator onto lead a and absorbed square root factors of the density of 
states in the leads into the coupling matrix W for notational simplicity. Note that we allow W to depend on energy. 
(Compare with the wide-band limit discussed in Ref. [22 , which employs an energy-independent hybridization F.) 
Dyson's equation for the retarded Green's function can then be written, in matrix form, as 

- ^9t'^^(^, t') = S{t - t') + / dtiG^'it, ti)I]«(ti, t') + g'^it, t')ho{X) . (17) 



To perform the adiabatic expansion, it is convenient to work in the Wigner representation, in which fast and slow 
time scales are easily identifiable. The Wigner transform of a function A(ti,t2) depending on two time arguments is 
given by 

A{t, e) = y dr e"^ A(t r/2, t - t/2) . (18) 

Using this prescription for the Green's function , the slow mechanical motion implies that varies slowly with 
the central time t = and oscillates fast with the relative time r = ti — t2- The Wigner transform of a convolution 

C{h,t2) - / dt3A{ti,t3)B{t3,t2) is given by 



C = exp 



AB 



~ is + '-d.AdtB - '-dtAd.B, (19) 

where we have dropped higher order derivatives in the last line, exploiting the slow variation with t. Therefore, using 
Eq. (19) we can rewrite the Dyson equation Eq. (17 1 as 

1 « (e - 1]« - ho) - '-d^g^'dtho ~ '-dtg"" (i - a.s^) , (20) 

where the Green's functions are now in the Wigner representation. Unless otherwise denoted by explicitly stating the 
variables, here and in the following all functions are in the Wigner representation. Finally, with the help of Eqs. ( A5 ) 
-(A6) from Supp. Mat. [Aj we obtain 

gi^c^G^'+'-Y, (9eG^A.G« - G^A,aeG«) , (21) 
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in terms of the strictly adiabatic Green's function 

G"ie,X)^[e~ho{X)-j:^(e)]- 



(22) 



Our notation is such that Q denotes full Green's functions, while G denotes the strictly adiabatic (or frozen) Green's 
functions that are evaluated for a fixed value of X (so that all derivatives with respect to central time in Eq. ( 20 ) 
can be dropped). From now on, (y(^i^><^>) denote the Green functions in the Wigner representation, with arguments 
(e,t),and^^ = (e«)t. 

Using Langreth's rule (see e.g. Ref. [15] ) 



(23) 



we can relate g< with g^. In Eq. (p3| we have introduced the lesser self energy S< , which in the Wigner representation 
takes the form 



E<(6)=2z^/„(6)r°(6) 



(24) 



Note that depends only on e and is independent of the central time. Expanding Eq. (231 up to the leading 
adiabatic correction according to Eq. (19 1, we obtain ^< to first order in X, 



g< = G< + - ^ i:^ [(9eG<)A^G^ - G"Kd,G< + (a,G^)A^G< - G<A^9,G^] , 
with G< = G^I]<G^. 



(25) 



C. Current-induced forces in terms of Green's functions 



We can now collect the results from the previous section and identify the current-induced forces appearing in the 
Langevin equation ([ij. Except for the stochastic noise force, the current induced forces are encoded in tr(tJ^A^). In 
the strictly adiabatic limit, i.e., retaining only the first term on the RHS of Eq. (251, g"^ ~ G^, we obtain the mean 
force 



de 



tr [A^G"^ 



(26) 



The leading order correction in Eq. ( 25 1 gives a velocity-dependent contribution to the current induced forces, which 



determines the tensor ^^u' ■ After integration by parts, we find 



J ^tr (G<A,aeG«A,, - G<A,,d,G^A,) 



This tensor can be split into symmetric and anti-symmetric contributions, 7 = 7^-1-7°, which define a dissipative 
term 7* and an orbital, effective magnetic field 7" in the space of the collective modes. The latter interpretation is 
based on the fact that the corresponding force takes a Lorentz-like form. Using Eq. (Al) in the Supp. Mat. [A| and 
noting that 2 / deG'^d^G'^ = J ded^{G'^)^ = 0, we obtain the explicit expressions 



7:.'(X) 
7.V(X) 

Here we have introduced the notation 



J gtr{A.G<A.,9,G>}^, 

-J gtr{A.G<A,a,(G^ + G^)}^ 



(27) 
(28) 



for symmetric and anti-symmetric parts of an arbitrary matrix A. 
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At last, the stochastic force is given by the thermal and non-equilibrium fluctuations of the force operator 
—SAiyd in Eq. ([t]). As indicated by the fluctuation-dissipation theorem, the fluctuating force is of the same order 
in the adiabatic expansion as the velocity dependent force. Thus, we can evaluate the expression for the correlator 
Dij^iCX.) of the fluctuating force given in Eq. (12) to lowest order in the adiabatic expansion, so that 



D,,,{X)= J ^tr{A,G<A,,G>}^ 



(29) 



This formalism gives the tools needed to describe the dynamics of the vibrational modes in the presence of a bias for 
an arbitrary number of modes and dot levels. When expressions (26) - ([28]) are inserted back in Eq. ([T]), they define 
a non-linear Langevin equation due to their non-trivial dependences on X(i) |28|29j . 



III. S-MATRIX THEORY OF CURRENT-INDUCED FORCES 

A. Adiabatic expansion of the S-matrix 

Scattering matrix approaches to mesoscopic transport generally involve expressions in terms of the elastic S-matrix. 
For our problem, the S-matrix is elastic only in the strictly adiabatic limit, in which it is evaluated for a fixed value 
of X, 

S{e,X) ^ I - 2TriW{e)G"{e,X)W'' (e) . (30) 



As pointed out by Moskalets and Biittiker |8l44j . this is not sufhcient for general out of equilibrium situations, even 
when X(t) varies in time adiabatically. In their work, they calculated, within a Floquet formalism, the leading 
correction to the strictly adiabatic S-matrix. We follow here the same approach, rephrased in terms of the Wigner 
representation. The full S-matrix can be written as [JS] (note that, in line with the notation established before for 
the Green's functions, the strictly adiabatic S-matrix is denoted by S while the full S-matrix is denoted by S) 

S{e, t) - 1 - 27ri [l^^;^M/t] (g^ (31) 
To go beyond the frozen approximation, we expand S to leading order in X, 

S{e, t) ~ S{e, X(<)) + ^ X,{t)A,{e, X{t)) . (32) 

V 

Thus, the leading correction defines the matrix A, which, similar to S', has definite symmetry properties. In particular, 
if the system is time-reversal invariant, the adiabatic S-matrix is even under time reversal while A is odd. For a given 
problem, the A-matrix has to be obtained along with S. 



We can now derive a Green's function expression for the matrix A |46I47| . Comparing Eq. ( 32 ) with the expansion 



to the same order of S in terms of adiabatic Green's functions (obtained straightforwardly by performing explicitly 

(33) 



the convolution in Eq. (31) and keeping terms up to X) we obtain 

A,(e,X) = ^9, [W^(e)G^(e,X)] A,(X)G«(e, X)iyt(g) 
-7rTy(e)G^(e,X)A,(X)9, [G«(e, X)W^t(g)] . 

Current conservation constrains both the frozen and full scattering matrices to be unitary. From the unitarity of the 
frozen S-matrix, S'^ S — 1, we obtain the useful relation 

^^'5 + 5t5^_0 (34) 



We will make use of Eq. (34) repeatedly in the following sections. On the other hand, unitarity of the full S-matrix, 



S'^ S — 1, imposes a relation between the A-matrix and the frozen S-matrix. To first order in the velocity X we have 

where A(e,X) — ^i^(e, X)Arjy. Therefore, S and A are related through 

A^S^ + SAl^- ——-——] . (36) 



2 \dX^ de de dX^ 

In the next section we will see that the A-matrix is essential to express the current-induced dissipation and "Lorentz 



forces, Eqs. (27) and (28) 
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B. Current-induced forces 



Mean Force 



The mean force exerted by the electrons on the oscihator is given by Eq. (26). Writing Eq. (26) exphcitly and 
using Eq. (A2) in Supp. Mat. [Aj we can express in terms of and and obtain 



i^.(X) = - / &^/„tr(A,G^Wtn„W^G^) 
= -J de^/„tr(n„W^G^A,G«M/t) , 



(37) 



where the second equahty exploits the cyclic invariance of the trace. Noting that, by Eq. (A7) in Supp. Mat. |Aj 



2ni dX^ 



Eq. (37) can be expressed directly in terms of scattering matrices S{e, X) as 

^''(X) = E/|^/"Tr(n„5t 



dS 



(38) 



(39) 



Note that now the trace (denoted by "Tr") is over lead-space. 

An important issue is whether this force is conservative^ i.e., derivable from a potential. A necessary condition for 
this is a vanishing "curl" of the force, 



dX^, 



E 



*f T /^n 



(40) 



From Eq. ( 40 ) it is seen that the mean force is conservative in thermal equilibrium, where Eq. ( 40 ) can be turned 



into a trace over a commutator of finite-dimensional matrices: Indeed, in equilibrium the sum over the lead indices 
can be directly performed since /« = / for all a, and Ila = 1. Using the unitarity of the S-matrix and the cyclic 
property of the trace, we obtain: 



f ds^ ds ds^ ^ ds 

2m \dX^dX^, dX^, dX^ 
— fT f ^ J^^\ -n 



(41) 



where in the last line we have used Eq. (34). In general, however, the mean force will be non- conservative in out-of- 



equilibrium situations, providing a way to exert work on the mechanical degrees of freedom by controlling the external 
bias potential [30148149) . 



2. Stochastic Force 



Next, we discuss the fluctuating force with variance D,^^' given by Eq. (29). Following a similar path as described 



in the previous subsection for the mean force F^, , we can also express the variance Eq. ( 29 ) of the fluctuating force 
in terms of the adiabatic S-matrix, 



de 



5t 



dS 
dX~, 



t 



dS 

ax; 



(42) 



where we have introduced the function F^a'i^) ~ /^(e) [1 — fa'{^)]- From Eq. (42) it is straightforward to show that 
D^^i is positive definite. By performing a unitary transformation to a basis in which D^^i is diagonal, using IIq, = 11^ 
and the cyclic invariance of the trace, we obtain the expression 



i?..(X) = 



OLOi' 



de 
2n' 



n<,.st|fn„ 

oX^ 



(43) 



which is evidently positive. 
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3. Damping Matrix 



So far, we were able to express quantities in terms of the frozen S-matrix only. This is no longer the case for the 
first correction to the strictly adiabatic approximation, given by Eqs. (27) and (28 1. We start here with the first of 



these terms, the symmetric matrix 7* , which is responsible for dissipation of the mechanical system into the electronic 
bath. 

The manipulations to write the dissipation term as a function of S-matrix quantities are lengthy and the details are 
given in the Supp. Mat. |B] The damping matrix can be split into an "equilibrium" contribution, ^'^•'^1^ and a purely 
non-equilibrium contribution 7'*'"'=, as 7" — 7'*''^'? + 7'*'"^. We first treat 7*^*^''. By the calculations given in Supp. Mat. 
IBI we obtain 



s.eq 



iE/|a</" + /-)Tv{n,s.J|n.5.^ 

aa' 



(44) 



where we have used that J2a'^a' = 1 , 5*^5 = 1, and Eq. (34) in the last line. Note that in general, 7^''^'? also 
contains non-equilibrium contributions, but gives the only contribution to the damping matrix when in equilibrium. 
Eq. ( 44 1 is analogous to the S-matrix expression obtained for dissipation in fcrromagnets in thermal equilibrium, 



dubbed Gilbert damping |35j . 

To express 7"'"*^ in terms of S-matrix quantities, we have to make use of the A-matrix defined in Eq. (33). Again 
the details are given in the Supp. Mat. [b1 where we find after lengthy manipulations that 



s,ne 



(45) 



This quantity vanishes in equilibrium, as can be shown using the properties of the S and A matrices. Since the sum 
over leads can be directly performed in equilibrium, expression (451 involves 



Tr 



dS'f X dS 



= -Tr 



-Tr 



dS 
dS 



SA], 



/ ds ds^ ^ds^ ds 







(46) 



2 ^' ydX^ \dX^. de de dX^ 

where we have used the unitarity of S and the cyclic invariance of the trace multiple times. In the first equality, we 



inserted S^S = 1 and used Eq. (34 1, the second equality follows by inserting the identity (36 1 and using again (34) 



Finally, combining all terms we obtain an S-matrix expression for the full damping matrix 7* 

95t dS 



7L,(X) = -^ J ^ae/„Tr|n„ 



dx, dx,, 



Al 



" dX, 



(47) 



Note that in equilibrium, by the relation —d^f — /(I — f)/T and using Eq. (34|, the fluctuating force D and 
damping 7* are related via 



as required by the fluctuation-dissipation theorem. 



(48) 



Following a similar set of steps as shown above for the variance Diy^i in Eq. (43), 7^J^^' has positive eigenvalues. On 
the other hand, the sign of 7^^T is not fixed, allowing the possibility of negative eigenvalues of 7^*. The possibility 
of negative damping is, therefore, a pure non-equilibrium effect. Several recent papers found negative damping in 
specific out of equilibrium models [22140150151] . 
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4- Lorentz force 



We turn now to the remaining term, the antisymmetric contribution 7" given in Eq. ( 28 ) , which acts as an effective 
magnetic field. Using Eq. (A2) in Supp. Mat. [A] it can be written as 



(49) 



In order to relate this to the scattering matrix, we use the Supp. Mat. |A]Eq. (AlOl, which allows us to write 7" in 
terms of the S-matrix as 



OA, 



(50) 



If the system is time-reversal invariant, 7° vanishes in thermal equilibrium. The latter implies ^^^afa = /, so 
that Eq. (50) involves only 



Tr 



9A]^ 5"* — ^^*^ 



dX,, 



Tr 



dA, . dAl 



yielding 7° = due to the cyclic invariance of the trace. In the last equality, we have used S ~ S"^ and A = —A^ as 
implied by time-reversal invariance. 

Out of equilibrium, 7° generally does not vanish even for time reversal symmetric conductors, since the current 
effectively breaks time reversal symmetry. 



IV. CURRENT 



So far we have focused on the effect of the electrons on the mechanical degrees of freedom. For a complete picture, we 
also need to consider the reverse effect of the mechanical vibrations on the electronic current. In the strictly adiabatic 
limit, this obviously has to reduce to the Landauer-Biittiker formula for the transport current. Considering the leading 
adiabatic correction to the current in equilibrium is closely related to the phenomenon of quantum pumping, and we 
will see that our results in this limit essentially reduce to Brouwer's S-matrix formula for the pumping current f^. 
Our full result is, however, more general since it gives the leading adiabatic correction to the current in arbitrary 
non- equilibrium situations (8|. 

The current through lead a is given by [43]: 

= -e(7V„) =ieY, W^n{cl{t)dn{t)) + h.c. (51) 

with Na — X^ijgQ c^c,,. Using the expressions for the self-energies this can be expressed in terms of the dot's Green's 
functions and self-energies, 

io.it) = e J dt' ti {g"{t,t')j:<{t' ,t) + g<it,t')j:i{t' ,t)} + h.c. . (52) 

Again we use the separation of time scales and go to the Wigner representation, yielding 

ia = ej ^tr{g^j:< + g<j:i-^{dtg^d,j:< + dtg<d,^i)} + h.c.. (53) 

We split the current into an adiabatic contribution and a term proportional to the velocity X^: 

Ic=ll + Il (54) 
We will express these quantities in terms of the scattering matrix. 
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Landauer-Biittiker current 



The strictly adiabatic contribution to the current is given by 

7° (X) = e I ^ tr { (G^ - G^^) S< + G< (E;^ - Sf ) } 



(55) 



where we have cohected the purely adiabatic terms from Eqs. (21) and (25). Inserting the expressions for the 
self-energies Eqs. (15) and (24), we can express this as 



(56) 



/°(X) = e / ^ ^/^2^*Tr{M/[5,;3(G«-G-4) + 27rzG«W^tn^M^G^]M^tn„}^ 

where we used Supp. Mat. |A]Eq. (A2). Inserting the adiabatic S-matrix, Eq. (30 1 yields 

/^(X) = e /■ ^^/,Tr{[<5„^-5n^5t]n„} (57) 

= e / ^ ^ (/„ - /^) Tr [SUpS^U^] , (58) 

^ /3 

where we used 511^5*^ = 1 in the last line. We hence recover the usual expression for the Landauer-Biittiker 
current [3]. Note that the total adiabatic current depends implicitly on time through X(i), and is conserved at every 
instant of time, ^aP^) — 0- To obtain the dc current, we need to average this expression over the Langevin 
dynamics of the mechanical degrees of freedom. Alternatively, we can average the current expression with the prob- 
ability distribution of X, which can be obtained from the corresponding Fokker-Planck equation. Similar remarks 
would apply to calculations of the current noise. 



B. First order correction 



Now we turn to the first order correction to the adiabatic approximation [5], restricting our considerations to the 
wide-band limit. The contribution to the current ( 53 ) which is linear in the velocity reads 



/i(X) = e J ^*^A:^tr{(a,G^)A^G^E< + [(a,G<) A^G^ - G^A^(a,G<)] } -t-h.c. 



(59) 



after integration by parts. Again, we insert Eq. (A2) from Supp. Mat. [Ajfor the lesser Green's function, and expres- 
sions ^ and ([24]) for the self-energies. In the wide band limit, the identity {i/2)d^dx,S + A,, = W{d^G^)A^G^W^ 
holds, so that we can write 



/;(x).-./|x^i:/,Tv[(l|j| + A)n,5.n„ 



(60) 



after straightforward calculation. After integration by parts, we can split this expression as 



/I (X) = -f / deX . ^ SJ^ImTr {^c^^pS^ 

deX ■ J2 U ReTr |jn„— n^j— - 2n„An;35^ 



(61) 



In equilibrium, the second term vanishes due to the identity Eq. (36) and the first term agrees with Brouwer's formula 
for the pumping current As for the strictly adiabatic contribution, the dc current is obtained by averaging over 
the probability distribution of X. 
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V. APPLICATIONS 



Resonant Level 



To connect with the existing Hterature, as a first example we treat the simplest case within our formalism: a resonant 
level coupled to a single vibrational mode and attached to two leads on the left (L) and right (R). This model has 
been discussed in detail for zero temperature in references |28I29| , and it provides a simple description on how current- 
induced forces can be used to manipulate a molecular switch. Here we derive finite-temperature expressions for the 
current-induced forces for a generic coupling between electronic and mechanical degrees of freedom, starting from 
the scattering matrix of the system, and show how they reduce to the known results for zero temperature and linear 
coupling. 

We consider N ^ M = 1, denoting the mode coordinate by X, the energy of the dot level by e(X), and the number 
of channels in the left and right leads by N]^ and Nfi, respectively. The Hamiltonian of the dot can then be written 
as 



Hr 



%x)d)d 



(62) 



and the hybridization matrix as W'^ = (w^, w^)^, with w" = (w", . . . w'^^) and a = L, R. Hence the frozen S-matrix, 
Eq. (30), is given by 



S 



2m 

IT 



(63) 



where £(e,X) = (-l{X)+iT, T = Tl+Tr, and = 7r(v^f")' -w". Rotating to an eigenbasis of the lead channels, this 
S-matrix does not mix channels within the same lead, and hence we can project the S-matrix into a single non-trivial 
channel in each lead, to obtain 



5 = 1 



2« / Tl ^/TlTr 



(64) 



To calculate the mean force from Eq. (39), we need an explicit expression for Eq. (A7) in Supp. Mat. |A] This can 
be easily calculated to be 



5t 



dS 
dX 



dX 1^ VTZTfl 



(65) 



and hence 



F{X) = 



de 

TT 



\c\' 



Analogously, the variance of the stochastic force, Eq. (42), becomes 



D{X) = 2 



dX 



dX 



(66) 



(67) 



It only remains to calculate the dissipation coefficient 7. Since there is only one collective mode, = 1, 7 is a scalar 
and hence 7° — 0. Moreover, for energy-independent hybridization we have that d^Gn — —G\, and the A-matrix 
( 33 ) can be written as [22] 



= -ttWGr[Gr,K]GbW'^ . 



(68) 



Being the commutator of scalars, in this case Ai — and from Eq. (47), 7^* must be positive and is given by Eq 



|44j). (For an alternative derivation of the positiveness of the friction coefhcient in a resonant-level system, see Ref. 
[52J). After some algebra, we obtain 



as 

dX 



dS 
dX 



= 4 



dX 



\C\ 



V^L^R ^R 



(69) 
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and hence the damping coefRcient becomes 



7(X) - - 



de IjAfL + TndjR 



de_ 
dX 



(70) 



We can evaluate the remaining integrals analytically in the zero-temperature limit |28l29j . In the following we 
assume > ^iR■ The average force is given by 



a 

Similarly we obtain the dissipation coefficient 



arctan 



7^(X) 



r 



dX 



E 



(71) 



(72) 



together with the fluctuation kernel 



D{X) 



7rr3 



dX 



arctan 



The position of the dot electronic level can be adjusted by an external gate voltage 



^J■L + ^J■R 



gate 



(73) 



(74) 



where the factor (/i^ + I^r)/'^ is included for convenience, to measure energies from the center of the conduction 
window. The difference in chemical potential between the leads is adjusted via a bias voltage 



eVbi; 



(75) 



For a single vibrational mode, the average current-induced force is necessarily conservative and we can define a 
corresponding potential. Restricting now our results to linear coupling, we write the local level as e{X) = cq + XX. 
In Fig. [ij we show the effective potential U{X) = ^UqX'^ — J dXF{X) which describes both the elastic and the 
current-induced forces at zero temperature and various bias voltages. Already this simple example shows that the 
current-induced forces can affect the mechanical motion qualitatively [29]. Indeed, the effective potential U{X) can 
become multistable even for a purely harmonic elastic force and depends sensitively on the applied bias voltage. 

Alternative expressions of the current-induced forces for the resonant level model, in terms of phase shifts and 
transmission coefficients, are given in the Supplementary Material [Cl 



B. Two-level model 



For the resonant level model discussed so far, the A-matrix vanishes and the damping is necessarily positive. We 
now consider a model which allows for negative damping jSSj. Our toy model could be inspired by a double dot on 
a suspended carbon nanotube, or an H2 molecule in a break junction. The model is depicted schematically in Fig. 
[2j The bare dot Hamiltonian corresponds to degenerate electronic states eo, localized on the left and right atoms or 
quantum dots, with tunnel coupling t in between, 

We consider a single oscillator mode with coordinate X that couples linearly to the difference in the occupation of the 
levels. In our previous notation, this means Ai — AiCs, where we denote by cr^, with /i = 0, . . . , 3, the Pauli matrices 
acting in the two-site basis. The shift of the electronic levels is given by e±{X) = eo ± \iX. The hybridization 
matrices are given by F" — ^Ta{cr^ ± cr^), where the -l-(— ) refers to a = L{R). We can deduce the tunneling matrix 
W in terms of the hybridization matrices, 



W = l/2yiV^((7" + G^) + l/2yiV^(a° - a^) 



(77) 
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FIG. 1: Resonant level. The shape of the effective potential U{X) can be tuned by the bias voltage. We consider the parameters 
eVgate = 0, hujo ~ 0.01 and F = 0.1. The dimensionless coordinate is a; = (AIluq/ X)X and energies are measured in units of 
AV(Ma;g). 




FIG. 2: Sketch of the two-level model. Electrons tunnel through two degenerate energy levels between left and right leads. The 
system is modulated by the coupling to the vibrational modes. 



In the wide-band limit, we approximate W and to be independent of energy. The retarded adiabatic GF takes the 
form 



(78) 



with A(X) = {e-i^+ iTL){e - e+ + iTR) - t^. 

For simplicity, we restrict our attention to symmetric couphngs to the leads, — — T/2. Hence the frozen 
S-matrix S{e,X) becomes 



while the A-matrix takes the form 



irfe-i++iT/2 tr 

A V e- +iT/2 



A{e,X) = iXiTt 



02 



(79) 



(80) 



We can now give explicit expressions for the current-induced forces. The explicit expressions are lengthy and are 
given in Supp. Mat. pi Eqs. (Dl| and (D2) for the mean force and damping matrix, respectively. The variance of 



the fluctuating force can be calculated accordingly. 

The average force given in Eq. (Dl) of Supp. Mat. combines with the elastic force to give rise to the effective 



potential U{X) depicted, for zero temperature, in Fig. k5l As in the case studied in the previous section, the system 
can exhibit various levels of multistability when changing the bias. 

The results for the friction coefficient, given in Supp. Mat. |D]Eq. (|D2| , are shown in Fig. |4] as a function of 
the dimensionless oscillator coordinate cc, for zero temperature. The contribution ^^^'^'i to the friction coefhcient is 
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FIG. 3: Effective potential for the mechanical motion in the two-level model. The shape of the potential can be tuned by 
changing the bias and gate voltages: (a) eVgato = 0, (b) elgatc ~ 0.2 and (c) eVgatc = 0.4. We consider the parameters 
hujo = 0.01, t = 0.1 and F = 0.1. The dimensionless coordinate is a; = {Mloq/\i)X and energies are measured in units of 
\1/{MloI). 



peaked at el^gato ± eVbias/2 = ±^ {XiXy + 1'^, as depicted in Figs, [i] (a) and (c). Neglecting the coupling to the 
leads, our toy model can be considered as a two-level system with level-spacing 2-y/ {XiXY + 1'^. Thus, the peaks 
occur when one of the dot's electronic levels enters the conduction window. When this happens, small changes in 
the oscillator coordinate X can have a large impact on the occupation of the levels. This effect is more pronounced 
when the dots' levels pass the Fermi levels that they are directly attached to [corresponding to X > for current 
flowing from left to right, see Fig. |4] (a) and Fig. [5] (a), (b)]. The broadening of the peaks is due to the hybridization 
with the leads, T/2. When eV^ate — 0, two peaks are expected symmetrically about X = 0, a,s shown in Fig. |4] 
(a) [see also Figs. Is] (a) and (b)]. The effect of a finite gate voltage eVgato is two-fold: it shifts the non-interacting 
electronic levels of the dot away from the middle of the conduction window, and hence the shifted levels e± pass the 
Fermi levels of right and left leads at different values of X, Figs. [5](c) and (d). Therefore in this case four peaks are 
expected, with two larger peaks located at X > 0, and two smaller peaks located at X < 0. This is shown in Fig. |4] 
(c). The height of the peaks in this case is reduced with respect to the case eV^ato = 0, since for a given peak, only 
one of the dot's levels is in resonance with one of the leads. Note that four real values of X can be obtained only if 
(eVgato ± eybias/2)^ > t'^. A situation with (eVgato - eVhias/^f < while (eVgato + eVbias/2)^ > is shown in 4[(c) 

(red-dotted line), where a big peak is observed for X — (eV^ato + eVbias/2)^ — t^, a corresponding small peak 

for X — —l/\i\J (eV^atc + eVbias/2)^ — t"^ [not displayed in Fig. |4j(c)], plus a peak dX X — Q. 

For this model, the A-matrix is generally non- vanishing, which can result in negative damping for out-of-equilibrium 
situations. This is due to a negative contribution of ^'^^^'^ to the total damping. This is visualized in Figs. [4] (b) 
and (d). Negative damping is possible when both dot levels are inside the conduction window, restricting the region 
in X over which negative damping can occur. Indeed, when only one level is within the conduction window, the 
system effectively reduces to the resonant level model for which, as we showed in the previous subsection, the friction 
coefficient is always positive. When current flows from left to right, negative damping occurs only for positive 
values of the oscillator coordinate X, as shown in Figs. |4](b) and (d). This is consistent with a level-inversion picture, 
as discussed recently in Ref. 51 . Pictorially, the electron- vibron coupling causes a splitting in energy of the left 
and right levels. When X > 0, electrons can go "down the ladder" formed by the energy levels by passing energy 
to the oscillator and hence amplifying the vibrations. For X < 0, electrons can pass between the two dots only by 
absorbing energy from the vibrations, causing additional non-equilibrium damping. For small broadening of the dot 



15 




FIG. 4: Damping vs. mechanical displacement in the two-level model, (a) Contribution ^"•'^'^ to the friction coefficient for 
various bias voltages at fixed gate voltage eV^atc = 0. (b) At the same gate voltage, the total damping exhibits a region of 
negative damping due to the contribution of 'y^-"'^. (c) 7"'^' for various gate voltages with the bias voltage eVbias ~ 0.8. Note 
that for both eV^ate = 0.2 and eVgato = 0.4, one small peak for negative x falls outside of the shown range of x. (d) Again, 
the full damping 7" exhibits regions of negative damping. We choose hcuo = 0.01, F = 0.1 and t = 0.1. The dimensionless 
coordinate is a:: = {Mu!q/\i)X and energies are measured in units of \i/{Mujq). 



levels due to the coupling to the leads, this effect is expected to be strongest when the vibration-induced splitting 
XiX becomes of the same order as the strength of the hopping t. When X grows further, the increasing detuning of 
the dot levels reduces the current and hence the non-equilibrium damping [see Figs. |4] (b) and (d) and Figs. [6] (a), 
(b)]. The coexistence of a multistable potential together with regions of negative damping can lead to interesting 
nonlinear behavior for the dynamics of the oscillator. In particular, and as we show in the next example, limit-cycle 
solutions are possible, in the spirit of a Van der Pol oscillator [M] ■ 

We can also calculate the current. The pumping contribution is proportional to the velocity X and thus small. 
Therefore we show here results only for the dominant adiabatic part of the current. This is given by 

For zero temperature, the behavior of the current is shown in Fig. 6] as a function of various parameters. Figs. [6] (a) 
and (b) show the current as a function of the (dimensionless) oscil ator coordinate x for two different values of gate 
potential for which the system exhibits multistability by developing several metastable equilibrium positions. For 
V^ate = and independently of bias, the current shows a maximum at the local minimum of the effective potential 
X = 0, while « for another possible local minimum, x w 0.5 (compare with Fig. |3](a)). The true equilibrium value 
of X can be tuned via the bias potential, showing the possibility of perfect switching. For finite gate potential however, 
the current is depleted from a; = with diminishing bias. Figs. [6] (c) to (d) show the current as a function of gate or 
bias voltage for fixed representative values of the oscillator coordinate x. The current changes stepwise as the number 
of levels inside the conduction window changes, coinciding with the peaks in the friction coefficient illustrated in Fig. 
[4j In an experimental setting, the measured dc current would involve an average over the probability distribution of 
the coordinate x, given by the solution of the Fokker-Planck equation associated to the Langevin equation ([T]). 
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(a) X > 



□ 



^gatc = 



(b) 



(d) 
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X <0 



■ Kate = 



a; < 



Vgate > 
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^^^^^ 



X > 



14ate > 



FIG. 5: Cartoon of the positions of the electronic levels in the dot with respect to the Fermi levels of the leads, depending on 
the sign of x and the existence of a gate voltage. The levels are broadened due to the hybridization with the leads F. When 
a; > 0, "left" and "right" levels approach the Fermi levels of left and right leads respectively, (a) for elgate = the levels align 
simultaneously for left and right, (c) a finite eVga.tc produces an assymmetry between left and right. For a; < the alignment 
of the levels is inverted, (b) eV^atc = 0, (d) finite eVgatc. 



C. Two vibrational modes 



As a final example, we present a simple model which allows for both a non-conservative force and an effective 
"Lorentz" force, in addition to negative damping. For this it is necessary to couple the two electronic orbitals of the 
previous example, see Eq. ( 76 1 , to at least two oscillatory modes which we assume to be degenerate. The relevant 
vibrations in this case can be thought of as a center-of-mass vibration Xi between the leads, and a stretching mode 
X2. (It should be noted that this is for visualization purposes only. In reality, for an H2 molecule, the stretching 
mode is a high energy mode when compared to a transverse and a rotational mode, see Ref. [55!. Nevertheless, the 
H2 molecule does indeed have two near-degenerate low energy vibrational modes, corresponding to rigid vibrations 
between the leads and a rigid rotation relative to the axis defined by the two leads.) The stretch mode modulates the 
hopping parameter, 

t ^ i{X2) = t + X2X2 , (82) 

while the center of mass mode Xi is modeled as coupling linearly to the density, 

eo^?(Xi) = eo + AiXi, (83) 

hence Ai = AiCTq and A2 = A2cri. We work in the wide-band limit, but allow for asymmetric coupling to the leads. 
The retarded Green's function becomes 

1 / e — e + iVji t 



where now A(A"i, X2) — {e — I + iT — e + i^R) — P- The frozen S-matrix can be easily calculated to be 

S{e,X,,X2)-l-^^ iVTZr^ ie~e + ^TI^)Tn) ■ ^^^> 
The A-matrices also take a simple form for this model. Since Ai is proportional to the identity operator, 

Aiie,Xi,X2) = -nXiWGR[GR,ao]GRW^ ^ 0. (86) 
On the other hand, the A-matrix associated with X2 is non-zero and given by 

A2(e,Xi,X2) = -iA2^^^(T2. (87) 
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FIG. 6: Dependence of the current in the two-level model on various parameters. Current as function of mechanical displacement 
for (a) Vgato = and (b) Vgatc = 0.4; as function of bias for (c) Vgato ~ 0, (d) Vgatc = 0.4, (e) a; = and (f) x — 0.5. We 
choose hujQ = 0.01, F = 0.1 and t = 0.1. The dimensionless coordinate is j: = {Mujfj\i)X and energies are measured in units 
of Xl/iMuji). 



From this we can compute the average force, damping, pseudo-Lorentz force, and noise terms. These are listed in 
Supp. Mat. [E] At zero temperature, it is possible to obtain analytical expressions for these current-induced forces. 
Studying the dynamics of the modes 2(i) implies solving the two coupled Langevin equations given by Eq. ([T]), 
after inserting the expressions for the forces given in Supp. Mat. [E] Within our formalism we are able to study 
the full non-linear dynamics of the problem, which brings out a plethora of new qualitative behavior. In particular, 
analyses which linearize the current-induced force about a static equilibrium point would predict run-away modes due 
to negative damping and non-conservative forces |30j . Taking into account nonlinearities allows one to find the new 
stable attractor of the motion. Indeed, we find that these linear instabilities typically result in dynamical equilibrium, 
namely limit-cycle dynamics |22) . We note in passing that limit cycle dynamics in a nanoelectromechanical system 
was also discussed recently in Ref. [53] . 

We have studied the zero-temperature dynamics of our two-level, two-mode system for different ranges of parameters. 
In Fig. ([t]) we map out the values of the curl of the mean force, (V x F)j_, indicating that the force is non-conservative 
throughout parameter space. We also plot one of the two eigenvalues of the dissipation matrix 7^, showing that it 
can take negative values in some regions of the parameter space. We find that it is possible to drive the system into 
a limit cycle by varying the bias potential. The existence of this limit cycle is shown in Fig. [s] (a), where we have 
plotted various Poincare sections of the non-linear system without fluctuations. The figure shows the trajectory in 
phase space of the (dimensionless) oscillator coordinate xi after the dynamical equilibrium is reached, for several cuts 
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FIG. 7: Curl of the average force and damping coefficient for the model with two vibrational modes: (a) The curl of the current- 
induced mean force F is, in a non-equilibrium situation, generally non-zero, indicating that the force is non-conservative, (b) 
One of the two eigenvalues of 7" . Remarkably, it undergoes sign changes. A dissipation matrix 7* which is non-positive definite 
implies destabilization of the static equilibrium solution found at lower bias potentials, in this case driving the system into a 
limit cycle, see main text and Fig. [s] The parameters used are such that A1/A2 = 3/2. The elastic modes are degenerate 
with fiojo = 0.014, Tl,r = ^^2 '^ {f^o ± o"^), and the hopping between the orbitals is t — 0.9. The dimensionless coordinates are 
Xi — {MujQ/X)Xi and energies are in units of /(Mloq), where A — (Ai + A2)/2. 




FIG. 8: Limit-cycle dynamics for the model with two vibrational modes, (a) At large bias voltages, Poincare sections of the 
four dimensional phase space show the presence of a limit cycle in the Langevin dynamics without fluctuating force, (b) Several 
periods of typical trajectories (for different initial conditions after a transient) in the presence of the fluctuating forces ^ are 
shown. The same general parameters as in Fig. [7] are used here. 



of the (dimensionless) coordinate X2- Each cut shows two points in xi phase space, indicating the entry and exit of 
the trajectory. Each point in the plot actually consists of several points that fall on top of each other, corresponding 
to every time the coordinate X2 has the value indicated in the legend of Fig. [S] (a). This shows the periodicity of the 
solution of the non-linear equations of motion for xi, X2 for the particular bias chosen. Surveying over the various 
values of X2 reveals a closed trajectory in the parametric coordinate space xi, X2- 

Remarkably, signatures of the limit cycle survive the inclusion of the Langevin force. Fig. [s] (b) depicts typical 
trajectories in the oscillator's coordinate space Xi, X2 in the presence of the stochastic force, showing fluctuating 
trajectories around the stable limit cycle. 

Experimentally, the signature of the limit cycle would be most directly reflected in the current-current correlation 
function, as depicted in Fig. [9) We flnd that in the absence of a limit cycle the system is dominated by two 
characteristic frequencies, shown by the peaks in Fig. |9] These frequencies correspond to the shift in energy of the 
two degenerate vibrational modes due to the average current-induced forces Fi and F2. When the bias voltage is 
such that the system enters a limit cycle, the current-current correlation shows instead only one peak as a function 
of frequency. This result, as shown in Fig. |9j is fairly robust to noise, making the onset of limit-cycle dynamics 
observable in experiment. 
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FIG. 9: Current-current correlation function in the presence of noise for the system with two vibrational modes. The limit cycle 
is signaled by a single peak {Vbias = 10, see Fig. [8|, as opposed to two peaks in the absence of a limit cycle (Vbiaa = 2.5, 5). 
Increasing the bias potential increases the noise levels but the peaks are still easily recognizable. The results are obtained by 
averaging over times long enough compared with the characteristic oscillation times. The same general parameters as in Fig. 
[3 are used here. 



VI. CONCLUSIONS 



Within a non-equilibrium Born-Oppenheimer approximation, the dynamics of a nanoelectromechanical system can 
be described in terms of a Langevin equation, in which the mechanical modes of the mesoscopic device are subject to 
current-induced forces. These forces include a mean force, which is independent of velocity and due to the average 
net force the electrons exert on the oscillator, a stochastic Langevin force which takes into account the thermal and 
non-equilibrium fluctuations with respect to the mean force value, and a force linear in the velocity of the modes. 
This last, velocity dependent force, consists of a dissipative term plus a term that can be interpreted as an effective 
"Lorentz" force, due to an effective magnetic field acting in the parameter space of the modes. 

In this work we have expressed these current-induced forces through the scattering matrix of the coherent mesoscopic 
conductor and its parametric derivatives, extending the results found previously in Ref. |22| . Our results are now valid 
for a generic coupling between the electrons and the vibrational degrees of freedom, given by a matrix /io(X), and 
for energy-dependent hybridization with the leads, given by the matrix W{e). We have shown that expressing all the 
current-induced forces in terms of the S-matrix is only possible by going beyond the strictly adiabatic approximation, 
and it is necessary to include the first order correction in the adiabatic expansion. This introduces a new fundamental 
quantity into the problem, the A-matrix, which needs to be calculated together with the frozen S-matrix for a given 
system. 

There are several circumstances in which the first non-adiabatic correction, encapsulated in the A-matrix, is nec- 
essary. While the average as well as the fluctuating force can be expressed solely in terms of the adiabatic S-matrix, 
the A-matrix enters both the frictional and the Lorcntz-like force. In equilibrium, the frictional force reduces to an 
expression in terms of the adiabatic S-matrix. Out of equilibrium, however, an important new contribution involving 
the A-matrix appears. In contrast, the A-matrix is always required to express the Lorentz-like force, even when the 
system is in thermal equilibrium. 

The expressions for the current-induced forces in terms of the scattering matrix allow us to extract important prop- 
erties from general symmetry arguments. Driving the nanoelectromechanical system out of equilibrium by imposing a 
bias results in qualitatively new features for the forces. We have shown that the mean force is non-conservative in this 
case, and that the dissipation coefficient acquires a non-equilibrium contribution that can be negative. We have also 
shown that when considering more than one mechanical degree of freedom, a pseudo Lorentz force is present even for 
a time-reversal invariant system, unless one also imposes thermal equilibrium on top of the time-reversal condition. 

Our model allows one to study, within a controlled approximation, the non-linear dynamics generated by the 
interplay between current and vibrational degrees of freedom, opening up the path for a systematic study of these 
devices. By means of simple model examples, we have shown that it is possible to drive a nanoelectromechanical 
system into interesting dynamically stable regimes such as a limit cycle, by varying the applied bias potential. In 
a limit cycle, the vibrational modes vary periodically in time, which can be the operating principle for a molecular 
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motor. On the other hand, the possibihty of non-conservative forces could also allow one to extract energy from the 
system, providing a controllable tool for cooling. The study of these kinds of phenomena for realistic systems is an 
interesting application of the formalism presented in this paper. 
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Appendix A: Useful relations 

Here we list a set of useful relations for the derivations in the main text. 



The Green's functions are related via 



1. Green's functions relations 



The lesser and larger Green's functions are given by 

Q. a. 

G> ^G< + G" - G-^ = -2TTi ^(1 - /„) G^M^^n^M^G-^. 



From ( 22 1 it is easy to see that 



and 



dx^G" = G^A^G^ 

d,G" = -G^(l - dj:")G". 
1. Green's functions and S-matrix relations 



Noting that (for given t) dx^G" = G^A^G^, we find using Eq.(A4): 



(9S 

S^^— = -27ri(l + 27riM^G-^M^^)W^G^A^G^W^^ = -27riW^G'^Aj.G^W^^ 



This holds for arbitrary magnitude of Xi, 
In the main text we use 



1 dS'f 
T^dxZ 



A^^ ^ 2mWG^A^G^W''d,{WG'^)A^'G"w'' - WG'^A^iG^ - G^)A^'9e(G^M^^), 



and 



5t 



dX,> 



= -2n [M^G'^A^(9eG^)A^,G^M^^]^ 



For energy- independent F" , we can use ( A6 ) so that also 



S^^^27nWG^G''w\ 
oe 



dS 



and ( A8 1 simplifies to 



dX, 



A,, = irWG^A, (G^ - G^) (A,,G^ - G^A,,) G^VF^ 



(Al) 

(A2) 
(A3) 

(A4) 
(A5) 

(A6) 



(A7) 

(A8) 
(A9) 

(AlO) 

(All) 
(A12) 

(A13) 
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Appendix B: S- matrix derivation of the damping matrix 



The expression for 7^* given in Eq. ( 27 1 can be written explicitly in terms of retarded and advanced Green's functions 
7^,, = 27r^y" deUY{A,G''W^UaWG^A,,d, [(1 - MG^'W^U^WG''] }^ . (Bl) 



We split Eq. (Bl) into two terms, the first due to the derivative acting on the Fermi function, the second from the 
rest, 7* = rys{i)~^js{ii) ^ r^Yjc flrst term is given by 



y^L'^ =2^E J deU~dJc,')tr{U^,WG''A,G''W^n^WG^A,,G''W^}^ 



(B2) 



where we have used the cyclic invariance of the trace. Similar to the derivation for the mean force, by means of 
expression (A7) in Supp. Mat. |Aj Eq. (B2) can be expressed in terms of the frozen S-matrix as 



The second contribution, in terms of G^- and G*^, reads 

= i^^)' E / {A^G^W^U^WG^A^,d, {G^W^U^,WG^)}^ 



(B3) 



(B4) 



It is instructive to split the factor F^a' into a symmetric and an antisymmetric part under exchange of the lead 
indices, F„„, = F^^, + F^^,, with 



(B5) 



Correspondingly, we split 7*^^^) into symmetric [7'*(^^'^)] and antisymmetric [7*(^^'')] parts in the lead indices: 7''(^-'^) = 
^s(iis) _|_ ^s{iia) ^ Due to its symmetries, 7'*^^^'*) can be easily expressed in terms of the S-matrix, 



s(IIs) 



= 7r^ j deF'^^,d,iT:{A^GVlicWG'^A^,G^W^naWG'^]^ 

QlOl' 

= -7r^ y" de (d,Fl^,)%T{A^G^W'^'^o.WG^A^,G^W'^Vi.^,WG^)^ 

OlOl' 



(B6) 



where in the second line we have integrated by parts since F^ vanishes for e — >■ ±00, and in the last line we have used 
Eq. (A7) from App. [A|once again. 



1. "Equilibrium" dissipative term 7" '^'' 
Since in equilibrium F'^^, — F^^^ = 0, 7*(^^°) I — and we can now regroup terms into an "equilibrium" contri- 

I eq 

bution, 7^''^9 = j^i^) + and a purely non-equilibrium contribution 7*^'"^ = ^s(^^i). 



(B7) 



By adding up expressions (B3) and (B6), it is straightforward to obtain Eq. (44) for 7'*''^' given in the main text. 
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2. Non-equilibrium dissipative term 7"'" 

To obtain 7^'"*^ in terms of S-matrix quantities we start from the expression 

s 

7^ 



and exploiting 11^ — 1 and the identity (A7) in Supp. Mat. [Aj we note that Eq. ( |B8[ ) can be written as 



s,ne 



de^/„trin„ 



^rQ 

dx,/ 



ae oe 



where [. , .] indicates the commutator. Calculating each term in the commutator separately we obtain 



(B8) 



(B9) 



5t 



dS 



de 



diG^W'f) d{WG' 



de 



A^-G'^W'' 



-WG^A^iG-^ ~ G")A^> 



djG^W'^) 
de 



de 



de 



-A.^G'^W^ 



+ 2^iWG^A,G''W^ d{WG^) ^^^^ 
de 

dS ^ JJWG^^^^^^^A _ G^)A,G«wt 
dXy de 

- 2^iWG^A,, ^(^^^^) WG^AM'^W^ 
de 



(BIO) 



where we have used Eq. (A4) from Supp. M at. [A} Finally, with help of the identity (A8) in Supp. Mat. [A] the 
non-equilibrium term can be expressed as Eq. (1451) in the main text. 



Appendix C: Resonant level forces: alternative expressions 



To calculate the current-induced forces for the resonant level model presented in Sec. |V| we can alternatively start 
with the popular S-matrix parametrization 



S 



(CI) 



where the transmission coefficient T and the phases 77, 9 depend on X. We present here the results for linear coupling, 
e{X) = eo + AX. We can then identify the transmission probability 



r{e,x) = 



{e-eo- xxy + r2 



(C2) 



and the phases 



r^ie,X) = -2--^tan(^— — ^ 
e{e,X) = |+, + arctan(-I|_I|_ 

We can now relate the current-induced forces to this S-matrix parametrization. The result for the average force can 
be split into a non-equilibrium force F"*^ and an equilibrium force F'^'^, i.e., F — F""^ + F'^'^ with 



(C3) 
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The amplitude of the fluctuating force can be obtained from Eq. ( 42 1 and is given by 



(C4) 



where we have defined 



Yll 



(1-r) 
(i-r) 





di] 


dX 


^ dX 






dX 


dX 


1 


fdT 



AT{l-T) \dX 



+ T{l-T) 



9(77 - 9) 



dX 



After some algebra, we also obtain 



2T 



(C5) 



This last expression corresponds to 7'*'^'? given in Eq. (44 1. (As we pointed out previously, ^"^'^'^ vanishes in this 
case). Here we have isolated a term that vanishes in equilibrium, showing explicitly that there is a non-equilibrium 
contribution in ( 44 ) . 



Appendix D: Current-induced forces for the two-level model 

The mean force is given by 



de 
2ti 



F{X) = -AiE 
The friction coefficient 7'' = r^^-'^i -\- ryS.ne j.gg^(jg 



2X,X{e~eo) {e - epf + jX^X)^ - + {T /2)' 

[JL + Jr) 1^|2 ^' [JL - JR) 1^|2 



(Dl) 



7 



s,eq 



An ■ — 



{2{e~eo)ty 



defR - OJl 



lAr 



((e - eo)2 + (r/2)2 + (AiA)2 + t^f + (2(e - eo)XiXf 
[4(e - eo)XiX ((e - eo)' + (r/2)2 + (X^X)^ - t')] \ , 



{ie-eo)'-{X,Xr-f)' 



+2(r/2)2 {{e - eof + (X^Xf + t') + {T/2)^ . 

Appendix E: Current-induced forces for the two vibrational modes model 



(D2) 



example discussed in the main text. For convenience, we define the following quantities: 

(6-6T + P + r?_ 



Here we list the current-induced forces quantities, calculated from Eqs. (|39|), (|42|), (|47|) and (50l for the two-modes 

(El) 
(E2) 



5ao(e) = 



lAI- 



2t(e- ?) 



5Q2(e) = ±- 



-2t Fi-Q, 



5a3(e) = ± — |2 



(E3) 
(E4) 
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where the +(-) refers to a = L{R) and with l-a = R{L) for a = L{R), and A(Xi,X2) = {e-e+iV L){e-e+iT r)-P ■ 



1. Mean force 



F. = 



-2 f ^AiV 

i 27r ^ [fg _ i)2 _ 12 _ 



/„(e)r„((e-e)2+P + rt„) 



~m2 



Fo. = - 



[{e - e)2 - P - TlTu] + [{Tl + rR)ie - e)] 
4 / t{^-~^){fL{e)rL + fR{e)TR) 



2. Fluctuating force 



(E5) 
(E6) 



Dii = 
Di2 = 
D22 = 



a0 



^ XI (1 - //3(e)) T/? (slaOS'/31 + 5al5/30) 

a/3 

2 (A2)^ J faie)^a (1 - //3(e)) (5a05/30 + 5al5;81 " 5c<25/32 - 9ai90i) 



(E7) 
(E8) 
(E9) 



3. Damping coefficients 



7ii 



7l2 = 



7I2 



(AiT 
27r 



/ de^(-5J„(e))rar0E 



/dc 



(ElO) 

(Ell) 
(E12) 



4. "Lorentz" term 



7l2 



/■ 



e — e 



(E13) 
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